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Abstract. In this paper we measure the two-body relaxation time from the angular deflection of test particles 
launched in a rigid configuration of field particles. We find that centrally concentrated configurations have re- 
laxation times that can be shorter than those of the corresponding homogeneous distributions by an order of 
magnitude or more. For homogeneous distributions we confirm that the relaxation time is proportional to the 
number of particles. On the other hand centrally concentrated configurations have a much shallower dependence, 
particularly for small values of the softening. The relaxation time increases with the inter-particle velocities and 
with softening. The latter dependence is not very strong, of the order of a factor of two when the softening 
is increased by an order of magnitude. Finally we show that relaxation times are the same on GRAPE-3 and 
GRAPE-4, dedicated computer boards with limited and high precision respectively. 
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1. Introduction 

CoUisionless iV-body simulations are heavily used for 
studying the dynamical evolution of galaxies, or systems 
of galaxies, and have so far been remarkably successful in 
producing many interesting results. Yet proper care has 
to be taken to eliminate possible sources of numerical er- 
rors which could, if present, lead to erroneous results. One 
of the possible sources of errors stems from the fact that 
the number of particles in a simulation is several orders 
of magnitude less than the number of stars in a typical 
galaxy, or, in other words, that the graininess in a com- 
puter realisation is much higher than that of the galactic 
system it is meant to represent. This could lead to errors 
since a particle moving through an A'^-body representa- 
tion of a given continuous system representing a galaxy 
is deflected from the orbit it would have had in the cor- 
responding smooth continuous medium, due to two-body 
encounters. This effect is known as two-body relaxation 
and the characteristic time linked to it as two-body relax- 
ation time (hereafter Trciax)- 

The relaxation time will obviously increase with the 
number of particles in the configuration, and tend to in- 
finity as the number of particles tends to infinity, in which 
limit the evolution of the system will follow the collision- 
less Boltzniann equation. Thus the relaxation time of sim- 
ulations will be much shorter than that of galaxies, which 
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is much longer than the age of the universe. Relaxation 
leads to a loss of memory of the initial conditions and an 
evolution of the system towards a state of higher entropy. 
It is thus necessary to have good estimates of the relax- 
ation times of iV-body simulations, since we can trust their 
results only for times considerably shorter than that. 

In the early times of A^-body simulations, when the 
number of particles used was of the order of a few hun- 
dred, the authors by necessity gave estimates of relaxation 
times in order to enhance the credibility of their results. 
Unfortunately in most cases only simple analytical esti- 
mates were used and the corresponding relaxation times 
were found to be comfortably, although perhaps unrealis- 
tically, high. As computers became faster, the number of 
particles used in simulations was increased. Authors using 
several tens or hundreds of thousands particles deemed 
unnecessary to include such simple estimates of the relax- 
ation times, since it was well known that the simple an- 
alytical estimates would give reassuringly high relaxation 
times. Nevertheless, it is not clear whether the simple an- 
alytical estimates are in all cases sufficiently near the true 
values. This could well be doubted since the simple analyt- 
ical estimates rely on a number of approximations, which 
are not in all cases valid. 

Since different A^-body methods may lead to different 
relaxation rates, it is of interest to discuss relaxation times 
when introducing a new method. It could t hus ha ve been 
feared that in a tree code (Barnes & Hut 1986 ) the re- 
laxation time would not be determined by the number of 
particles, but by the number of nodes, which would then 
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act as "super-particles" . Since the number of nodes is al- 
ways much smaller than the number of particles this would 
entail considerably shorter relaxation times than direct 
summation with the same number of particles, and thus 
constitute a major disadvantage of the tree code. This fear 
was put to rest by Hernquist (1987) who showed that the 
relaxation for tree code calculations does not differ greatly 
from that obtained by direct summation provided the tol- 
erance parameter is less than 1.2. Similarly Hernquist & 
Barnes ( p^90| ) compare relaxation rates in direct summa- 
tion, tree and s pherical harmonic A^-body codes, while 
Weinberg (|1996| ) introduces a modification of the orthog- 
onal function potential solver that minimises relaxation. 

Several methods have been us ed to measure two-body 
relaxation. Standish & Aksnes ( 1969 ), Lecar & Cruz- 



Conzalez ( 1972 ) and Hernquist ( 1987 ) have measured the 
angular deflection of test particles moving in a configura- 
tion of N field particles. Although this method has the 
disadvantage of not including collective effects, it has the 
advantage that all the parameters can be changed inde- 
pendently of each other, and that the results are easy to 
interpret. Theis (1998) performed semi-analytical calcula- 
tions, assuming a homogeneous medium and also ignor- 
ing cumulative effects. The most widely used approach is 
to monitor the energy conservation of individual particles 
in systems in which, had it not been for the individual 
encounters, the individual energies would have been con- 
served. This method, which includes collective effects, is 
well suited for testing relaxation rates introduced by dif- 
ferent codes, but can only be used with systems in equilib- 
rium. It has been used e.g. by Hernquist & Barnes ( 1990D , 



Hernquist & Ostriker ( 1992 ) , Huang, Dubin ski & Carlberg 
(|1993|) and Weinberg (|l996|) . Theuns (|l996|) measured the 



diffusion coefhcients as a function of the energy in a direct 
summation A^-body simulation by studying the properties 
of the random walk in energy space for particles of given 
energy and found very good agreement with theoretically 
calculated diffusion coefficients. Finally a number of stud- 
ies (e.g. Farouki & Salpeter |1982| , |1994| ; Smith |1992D rely 
on a measurement of the mass segregation, i.e. on the fact 
that, due to two-body encounters, high mass particles lose 
energy and spiral towards the center, while light ones gain 
energy and move to larger radii. Thus the configuration of 
the high mass particles contracts, while that of the hght 
particles expands, and from the rate at which this happens 
we can calculate the two-body relaxation time. 

In this paper we will calculate the relaxation times 
in a large number of cases, using the first of the meth- 
ods mentioned above, i.e. by measuring deflection angles 
of individual trajectories of test particles in a configura- 
tion of rigid field particles. We will cover a much larger 
part of the parameter space than was done so far, and 
we will also extend to larger number of particles. All the 
calculations presented in this paper were made on the 
Marseille GRAPE-3, GRAPE-4 and GRAPE-5 systems. 
The Marseille GRA PE-3 s ystems have been described by 
Athanassoula et al. ( 199 j ), while a general description of 
the GRAPE-4 systems and their PCI interfaces has been 



given by Makino et al. (1997) and Kawai et al. (1997) re- 
spectively. A description of the GRAPE-5 board can be 
found in Kawai et al. (2000). Opting for a GRAPE sys- 
tem restricts us to a single type of softening, the standard 
Plummer softening, but has the big advantage of allowing 
us to make a very large number of trials, covering well the 
relatively large parameter space. Theis (1998) compared 
the relaxation rates obtained with the standard Plummer 
softening to those given by a spline (Hernquist & Katz 



1989) and showed that the differences between the two 



are only of the order of 20-40%. 

This paper is organised as follows: In section ^ we 
briefly summarise the simple analytical estimates of the 
relaxation time. In section ^ we describe the numerical 
methods used in this paper and discuss the validity of 
their approximations. Here we also introduce the mass 
models which will be used throughout this paper. The 
values of the parameters to be used, and in particular the 
values of the softening, are derived and discussed in sec- 
tion ^. In section ^ we give results for the relaxation time. 
We specifically discuss the effect of number of particles, 
of the velocity and of the softening, and compare results 
obtained with GRAPE-3 and GRAPE-4. We also give a 
prediction for the relaxation time in an A'^-body simula- 
tion. We summarise in section ^ 

2. Simple analytical estimates of the relaxation 
time 

The relaxation time for a single star can be defined as 
the time necessary for two body encounters to change its 
velocity, or energy, by an amount of the same order as 
the initial velocity, or energy, i.e. the time in which the 
memory of the initial values is lost. Thus for the velocity 
we have 



Tr 



relax 



'Al 



(1) 



where Tcross is the crossin g time of the system. Following 
e.g. Binney & Tremaine (1987) we can obtain a simple 
order-of-magnitude estimate of the relaxation time. Let 
us focus on the motion of a single star assuming that it 
is moving on a straight line with a constant velocity v 
and that the remaining stars have equal mass m and are 
distributed uniformly in space. We first consider that the 
star passes a single perturber star on its rectilinear tra- 
jectory. Then the total change of its velocity component 
perpendicular to the trajectory is 

\Sv^\-^^, (2) 

where G is the gravitational constant and b is the impact 
parameter, i.e. the minimum distance between the two 
stars if there was no gravitational attraction. To find the 
total change, due to all the particles, we integrate over all 
encounters and find 



ln(T ), 



(3) 
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where N is the number of stars, m is their mass, R is 
some characteristic radius of the system and bmin and 5max 
are the minimum and maximum values for the impact 
parameter. Using the same approximations and including 
a softening e, Huang, Dubinski & Carlberg (1993) find 



R^v^ 



&2 



^min ' ^ 



62 
In "'^^ 



(4) 



A somewhat more elaborate derivation foUow- 
mg the calculation of the diffusion coefficients (e.; 



Chandrasekhar |1942|, Spi tzer & Hart |1971| , Spitzer |1987 
Binney & Tremaine 1987) gives 



^rclax — 



G^mp ln(fe„iax/6min) ' 



(5) 



where p is the density and _F is a constant, roughly equal 
to 0.34. The quantity ln(&max/6min) is often denoted by 
In A and referred to as the Coulomb logarithm. 

The appropriate values of fomin and femax in the above 
equations have b een d iscussed at length in the literature. 
Chandrasekhar (1942) argued that 6„iin is the value of b 
for which the angular deflection of the star is equal to 
iyr. The value of femax has been subject to considerable 



controversy. Chandrasekhar (1942), Kandrup (1980) and 
Smith (1992) have opted for a fomax of the order of the 
mean inter-particle distance, while others (e.g. Spitzer & 
Hart |1971| , Farouki & Salpeter |1982| , Spitzer |1987D used for 
bmax a characteristic radius of the system. The numerical 



simulations of Farouki & Salpeter (1994) argue in favour 
of the latter. This is further corroborated by the results 



of Theis (1998) 



Using the estimates 5niin = 6max = R and assum- 
ing virial equilibrium, so that we can use for the velocity 



the estimate 



19871 ) 
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we find (Binney and Tremaine 



(6) 



Similarly for the case with softening setting 6i„in ^ 
bmax = R and estimating the velocity by assuming virial 
equilibrium we find 



v^R^ 

G2M2 4[ln(i?2/e2) _ i] 

u4i?2 N 
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relax 
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T 
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-Tr, 



(7) 



Equation (Q) differs by a fact or of 2 from that of 
Huang, Dubinski & Carlberg (1993), the reason being that 
our definition of the relaxation time is based on t he vel oc- 
ities, while that of Huang, Dubinski & Carlberg (1993) is 
based on the energies. 



It is clear that, for the number of particles used in 
present day simulations of coUisionless systems and the ap- 
propriate values of the softening, N is considerably larger 
than R/e. Since, however, only the logarithms of these 
quantities enter in equations (|^) and (|7|) , the differences in 
the estimates of the relaxation times differ, for commonly 
used values of N, by less than a factor of 2. Equation ^ 
is more appropriate, since it includes the softening. Often 
a coefficient g is introduced in the Coulomb logarithm, 
i.e. A = gN. Giersz & Heggie (1994) estimated that the 
most appropriate value of g is 0.11. They also compiled 
in their Tabic 2 the values given by several other authors. 
They are all between 0.11 and 0.4. Independent of what 
is chosen for the Coulomb logarithm, equations such as 
(^ or (0) argue that even for a moderately low number of 
particles, of the order of say a few thousands, the relax- 
ation time is comfortably high, of the order of, or higher 
than, 40 crossing times. 

3. Numerical methods 

3.1. Method 



Following Standish & Aksnes ( 1969|) , Lecar & Cruz- 
Conzalez ( 1972 ) and Hernquist ( 1987 ) we will measure 
relaxation using the angular deflection suffered by a test 
particle moving in a configuration of N field particles fixed 
in space. 

Let us consider a sphere of a given density profile rep- 
resented by iV fixed particles, which we will hereafter refer 
to as field particles, and let us place a test particle of zero 
mass at the edge of this sphere either at rest, or with a 
radial velocity v. In the limit of N —^ oo its trajectory will 
be a straight line passing through the center of the sphere, 
which we will hereafter call the theoretical trajectory. For 
a finite N, however, the test particle will be deflected by 
a number of encounters with the field particles and thus 
it will cross the surface of the sphere at an angle $ from 
the corresponding theoretical point. Following Standish & 



Aksnes (1969), or Lecar & Cruz-Conzalez (1972), we can 
measure the relaxation time as 



'^rclax — 



Tt 



sm2<I> ' 



(8) 



where Tt is the crossing or transit time. Different realisa- 
tions of the adopted model will of course give somewhat 
different relaxation times. It is thus better to use an esti- 
mate based on the average of many realis ations or trajec- 
tories. Thus we have (Standish & Aksnes 1969 ) 



^rclax — 



<Tt> 
< sm2$ > ' 



(9) 



where the <> denote an average over all realisations 
and/or all trajectories. 

We will repeat such calculations here, extending them 
to non-homogeneous density distributions, different initial 
velocities of the test particle and a larger range of field 
particle numbers N. This will allow us to discuss the effect 
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of central concentration, of initial test particle velocities, 
of softening and of particle number on the relaxation time. 

3.2. Mass distributions and computing miscelanea 

To find the effect of central concentration on the re- 
laxation time we use three different mass distributions, 
namely the homogeneous sphere (hereafter model H), the 
Plummer model (here after model P) and the Dehnen 



sphere (Dehnen 1993 - hereafter model D). For the 



Plummer model the density is 
and for the Dehnen one 



p(r) 



(3 - 7)Mt 



47r 



r-f {r + a^Y'^-'ry' 



(10) 



(11) 



where Mt is the total mass of the sphere, ap and od are 
the two scale-lengths, of the Plummer and Dehnen models 
respectively, and 7 is the concentration index of the latter. 
For the model P we have taken a scale-length of ap = 9.2, 
and for model D we have adopted ap = 3.1 and 7 = 1. In 
all the examples in sections 5.1 to ^.4| we have truncated 



the mass distributions, using a cut-off radius of 20, and 
taken the mass within this cut-off radius to be equal to 
15. For these parameters the mass inside the cut-off radius 
is 75% of the total. For the homogeneous sphere we have 
adopted a cut-off radius of 20 and a mass of 15, so as to 
keep as much in line as possible with the other two models. 
For the gravitational constant we use G — 1. 

These three models span a large range of central con- 
centrations. They have the same cut-off radius, but the 
radius containing 1/10 of the total mass is for model D 
roughly one fourth of the corresponding radius for model 
P and roughly an order of magnitude smaller than that 
of model H. Similarly the radius containing half the mass 
for model D is roughly half of that of model P and a third 
of that of model H. These models will thus allow us to 
explore fully the effect of central concentration on the re- 
laxation time. 

We start 1 000 test particles from random positions on 
the surface of a sphere of radius 20 and with initial radial 
velocities equal to /v times their theoretical escape veloc- 
ity (hereafter Vcsc), calculated from the model. The parti- 
cles were advanced using a leap-frog scheme with a fixed 
time-step of At = 2~^ = 0.015625. We made sure this 
gives a sufficient accuracy by calculating orbits without 
the use of GRAPE and with this time-step. This showed 
that the energy is conserved to 8 or 9 digits. We then re- 
peate d the exercise using a Bulirsch-Stoer scheme (Press 
et al. 1992| ) and found that the energy was conserved to 10 
digits. Since an accuracy of 8 or 9 digits is ample for our 
work, we adopted the leap-frog integrator and the afore- 
mentioned time-step. 

The forces between particles were calculated using one 
of the Marseille GRAPE-3AF systems, except for the re- 
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Fig. 1. Distribution of the number of test particles orbits 
that have a given deflection angle, n($), as a function 
of that angle for model D and /v — 1.5, e = 0.01 and 
TV = 4 000 (a), iV = 2 000 (b) or = 1000 (c). 



suits given in sections 5.4 and p.5| , where we used the 
Marseille GRAPE-4 system. 

For each model we take 10 different field particles dis- 
tributions and for each we calculate the 1 000 test particles 
trajectories. For simplicity the test particles are the same 
in each of the 10 field particles distributions. For each of 
the test particles we calculate Tt and the deflection an- 
gle from the theoretical (straight line) orbit, <&. Then the 
relaxation time is calculated using equations (^) and (^. 
The errors are obtained using the bootstrap method (Press 
et al. 1992). In Figures || to || error bars are plotted only 
when aT^.i^^/Troiax > 0.05. 

3.3. Validity of the approximations 

Equations (^) and (^ were derived assuming small deflec- 
tion angles. There could, however, be cases in which a test 
particle comes very near a given field particle and is very 
strongly deviated from its initial trajectory, so that the de- 
flection angle is greater than 90°. In that case equations 
(H) and (H) are certainly not valid, particularly since for 
a deflection of 180° they will give the same result as for 
0°. It is not easy to treat such deflections accurately, so 
what we will do here is to keep track of their number and 
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make sure that it is sufficiently small so as not to influence 
our results. We have thus monitored the number of orbits 
whose velocity component along the axis which includes 
the initial radial velocity changes sign. We will hereafter 
for brevity call such orbits looping orbits. None were found 
for the homogeneous sphere distribution, and very few, 3 
in 10 000 at the most, for the case of the Plummer sphere, 
and that for the smallest of the softenings used here (cf. 
section 5.3 and Figure H). The largest number of loop- 
ing orbits was found, as expected, for model D and the 
smallest softening, i.e. e ~ 2^^^. For this case, /v = 1.5 
and N — 64000, we find of the order of 30 such orbits in 
10 000. Although this is considerably larger than the corre- 
sponding number for homogeneous and Plummer spheres, 
it is still low enough not to influence much our statistics, 
particularly if we take into account that it refers to an ex- 
ceedingly small softening. For the more reasonable value 
e = 2~^, we find that there are no looping orbits at all. 
We can thus safely conclude that the number of particles 
with looping orbits is too low to influence our results. 

We still have to make sure that the remaining orbits 
have sufficiently small deflection angles for equations (||) 
and ^ to be valid. Figure |l| shows a histogram of the 
number of test particles orbits that have a given deflection 
angle, n(<i>), as a function of that angle, <&, for model D, 
i.e. the one that should have the largest deflection angles, 
with /v = 1.5, e = 0.01. The upper panel corresponds to 

= 4 000, the middle one to A^ = 2 000 and the lower 
one to A^ = 1 000. Note that we have used a logarithmic 
scale for the ordinate, because otherwise the plot would 
show in all cases only a few bins near 0°. 

For A^ = 1 000 there is one particle with deviation of 
28°, one with 26° and two with 25°. Furthermore only 
30 trajectories, of a total of 10 000, have a deflection angle 
larger than 20° . The numbers are even more comforting for 
A^ — 2 000, where only two trajectories have a deflection 
angle larger than 20°, and even more so for A^ = 4 000, 
where no particles reach that deflection angle. We can thus 
conclude that in all but very few cases the small deflection 
angle hypothesis leading to equations (||) and (||) is valid. 



4. Free parameters 

We will calculate the values of the relaxation time as a 
function of three free parameters, namely the number of 
field particles, the initial velocity of the test particles and 
the softening. In this section we will discuss what the rel- 
evant ranges for these parameters are. 



4.1. Number of particles 

We will consider numbers of particles between 1 000 and 
64000. Indeed fewer than 1000 particles are hardly used 
anymore in numerical simulations, even with direct sum- 
mation. For more than 64 000 particles the two-body re- 
laxation is small. Furthermore the range considered is suf- 
ficiently large for all trends to be clearly seen. 



4.2. Initial velocity 

The question for the initial velocity is more convoluted. 
The simple analytical approaches leading to equations ^ 
and (0), instead of taking a spectrum of velocities for the 
individual encounters and then integrating over these ve- 
locities, introduce an effective or average velocity and as- 
sume that all interactions are made at this velocity. In 
our numerical calculations individual encounters occur at 
different velocities, depending on their position on the tra- 
jectory of the test particle and on the initial velocity of 
this particle. We can, nevertheless, introduce Vcs, an aver- 
age or effective velocity, in a similar way as the analytical 
approximation. A simple and straightforward, albeit not 
unique, such definition can be obtained as follows: Let us 
assume that the test particle moves on a straight line. 
At each point of its trajectory we can define a thin sheet 
going through this point and locally perpendicular to the 
trajectory. It will contain all field particles which have this 
point as their closest approach with the test particle. Let 
dr be the thickness of this sheet and Xdr the fraction of 
the total mass of the field particles that is in it. Then we 
can define the effective velocity 



VcS 



Xv dr. 



(12) 



which of course depends on both the distribution of the 
field particles and the initial velocity of the test particle. 

4.3. Softening 

The third free parameter in our calculations is the soft- 
ening. Merritt (199 6; herea fter M96 ) and Athanassoula et 
al. (2000; hereafter AFLB ) showed that, for a given mass 
distribution and a given number of particles A'^, there is a 
value of the softening, called optimal softening egpt, which 
gives the most accurate representation of the gravitational 
forces within the A^-body representation of the mass dis- 
tribution. For values of the softening smaller than Copt the 
error in the force calculation is mainly due to noise, be- 
cause of the graininess of the configuration. For values of 
the softening larger than Cgpt the error is mainly due to 
the biasing, since the force is heavily softened and there- 
fore far from Newtonian. Since the two-body relaxation is 
also a result of graininess it makes sense to consider soft- 
ening values for which it is the noise and not the bias that 
dominates, i.e. to concentrate our calculations mainly on 
values of the softening which are smaller than or of the 
order of Copt, keeping in mind that too small values of 
the softening have no practical significance. The value of 
the optimal softening decreases with A^ and can be well 
approximated by a power law 



^opt 



BN" 



( |M96| ) . The values of B and b depend on the mass dis- 
tribution under consideration, and, to a mu ch smaller de- 
gree, on the range of A'^ considered (lAFLBl) . Thus denser 
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configurations require smaller softenings for an optimum 
representation of the force (AFLB). For a given number 
of particles the homogeneous and Plummer spheres have 
roughly the same optimal softening, while the Dehnen 
sphere with 7 = has an optimal softening 0.45 dex lower 
(cf. Figure 9 ofjAFLB|). 

M96 and AFLB have calculated eopt using the mean 
average square error (MASE) of the force, which mea- 
sures how well the forces in an A^-body representation of 
a given mass distribution represent the true forces in this 
distribution. The average square error (ASE) is defined as 



ASE 



C_ 

N 



N 



(13) 





Fig. 2. MASE as a function of the softening e for the 
three models described in this paper. The upper panel 
gives the results for model H, the middle one for model P, 
and the lower one for model D. In each panel from top to 
bottom the curves correspond to = 1000, 2 000, 4 000, 
8 000, 16 000, 32 000 and 64 000, where N is the number of 
particles in the realisation of each model. The number of 
realisations taken in all cases is 6 xlO^/N. The position of 
the minimum error along a line corresponding to a given 
A'^ is marked by a x , and the corresponding e value is the 
optimal softening eopt for this number of particles. 



where Ftruo(xi) is the true force from a given mass distri- 
bution at a point Xj , is the force calculated at the same 
position from a given A^-body realisation of the mass dis- 
tribution and using a given softening and method, A^ is the 
number of particles in the realisation, and the summation 
is carried out over all the particles. In order to get rid of 
the dependence on the particular configuration, which is 
of no physical significance, many realisations of the same 
smooth model must be generated and the results should 
be averaged over them. Thus MASE, the mean value of 
the ASE, is 



r ^ 

MASE = ^ < E 

1=1 



-Ftruc(Xj)P > 



(14) 



where <> indicates an average over many realisations. In 
equations (|l^) and (|l^) C is a multiplicative constant, in- 
troduced to permit comparisons between different mass 
distributions. Since in this paper we are only interested in 
the values of eopt we will simply use C = 1. The MASE val- 
ues were found using 6 x 10^/A^ realisations and were cal- 
culated using direct summation on the Marseille GRAPE- 
5 systems. 

Figure ^ shows values of MASE as a function of e for 
the three models considered in subsections 5.1 to 5.4 and 



seven values oi N , in the range of values considered here 
(cf. subsection 4.1). The general form of the curves is as 
expected. There is in all cases a minimum error between 
the region dominated by the noise (small values of the 
softening) and the region dominated by the bias (large 
values of the softening). This minimum - marked with a 
X in Figure ^ - gives the value of eopt ■ For all three models 
a larger number of particles corresponds to a smaller error 
and a smaller value of Eopt, as expected (M96, AFLB). 

The more concentrated configurations give smaller val- 
ues of Eopt, as already discussed in AFLB, Thus for A^ = 



64000 the optimum softening for model H is less than 
twice that of model P, while the ratio between the soften- 
ings of models P and D is more than 10^ 



Comparing our results to those of [AFLBj we can get 
insight on the effect of the truncation radius. For this it 
is best to compare our res ults obtained with A^ = 32 000 
with those given by AFLB for A^ = 30 000, since these two 
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N values are very close and we do not have to make correc- 
tions for particle number. For our model P and = 32 000 
we find an optimum softening of 0.52, or, equivalently, 
O.OSTap, where ap is the scale length of the Plummcr 
sphere. This is smaller than the value of 0.063ap found for 
the lAFLBl Plummer model and the difference is due to the 

trun- 



AFLB 



different truncation radii of the two models, 
cated their Plummer sphere at a radius of 38.7fap, which 
encloses 0.999 of the total mass of the untruncated sphere, 
while model P is truncated at 2.2ap, which contains only 
75% of total mass. The difference in the values of Copt is 
in good a greement with the discussion in subsection 5.2 
of AFLB , When the truncation radius is large, the model 
includes a relatively high fraction of low density regions, 
where the inter-particle distances are large. This is not the 
case if the truncation radius is much smaller, as it is here. 
Thus the mean inter-particle distance is larger in the for- 
mer case and, as can be seen from Fig. 11 of AFLB , the 
corresponding optimal softening also. This predicts that 
the optimal softening should be smaller in models with 
smaller truncation radius, and it is indeed what we find 
here. 

Our model H is the same as that of AFLB , except 
for the difference in the cut-off radii. Thus the values of 
the optimum softening are the same, after the appropri- 
ate rescaling with the cut-off radii has been applied. Our 
model D differs in two ways from the corresponding model 
of |AFLE| . We use here 7=1, while |AFLE| used 7 = 0. 



We also truncated our model at 6.5aD, while AFLB trun- 
cated theirs at 2998aD- It is thus not possible to make any 
qualitative or quantitative comparisons. 



5. Relaxation measured from angular deflections 



5.1. The effect of the number of particles and of the 
density distribution 

Figure H shows the relaxation time as a function of the 
number of particles for the models H, P and D described 
in the previous section. The upper panel was obtained 
for /v = 1.5 and e = 0.01 and the lower one for /v = 
1.5 and e = 0.5. As can be seen from Figure ||, for the 
first value of the softening noise dominates over bias for 
all three models. For the second value noise dominates for 
model H, bias dominates for model D, and we are near the 
optimum value for model P and the high N values. The left 
ordinate in fig. ^ gives the relaxation time while the right 
one gives the ratio of the relaxation to transit time for 
the homogeneous sphere and /v — 1.5. The corresponding 
values of this ratio for the other two density distributions 
can be easily obtained if one takes into account that for 
/v = 1.5 the three theoretical transit times are 20.35, 17.52 
and 16.91, for the H, P and D models respectively. We 
fitted straight lines 

logio(T,.elax) ^Ai+Bi logio(7V) 



Table 1. Coefficients of linear regression in log-log scale 
of T,. 

oiax as a function of particle number N 



Model 


/v 


e 




Bi 


0"Ai 




H 






0.78 


1.00 


0.05 


0.01 


P 


1.5 


0.01 


0.76 


0.98 


0.04 


0.01 


D 






0.63 


0.78 


0.04 


0.01 


H 






0.87 


1.01 


0.05 


0.01 


P 


1.5 


0.03 


0.80 


1.00 


0.11 


0.03 


D 






0.86 


0.77 


0.11 


0.03 


H 






0.99 


1.01 


0.06 


0.01 


P 


1.5 


0.10 


0.94 


1.00 


0.11 


0.03 


D 






0.75 


0.89 


0.07 


0.02 


H 






1.13 


1.02 


0.08 


0.02 


P 


1.5 


0.50 


1.24 


0.98 


0.07 


0.02 


D 






1.29 


0.94 


0.07 


0.02 


H 






1.17 


1.01 


0.04 


0.01 


P 


2.0 


0.03 


0.98 


1.01 


0.08 


0.02 


D 






0.78 


0.87 


0.12 


0.03 


P 


3.0 


0.01 


1.40 


0.97 


0.04 


0.01 





o 
o 
o 



1000 



10" 



o 



o 



Fig. 3. Relaxation time as a function of the number of 
field particles, for /v=1.5 and three different mass mod- 
els, namely model H (stars), model P (x) and model D 
(crosses). The upper panel (a) was obtained with e = 0.01 
and the lower (b) with e — 0.5. The straight solid lines 
are the corresponding linear least square fits. The dashed 
lines give the predictions of equation (^, when we use 
^max — Ri while the dotted lines give the prediction of the 
corresponding equation obtained by using femax — I- 
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to each model using least-square fits and plotted them 
with solid lines in the figure. The values of Ai and Bi, as 
well as of their corresponding uncertainties aAi and ctbi, 
are given m Table |l|. We have performed similar calcula- 
tions for other values of the velocity and softening param- 
eters and have included in this Table the corresponding 
values of Ai and _Bi, as well as their uncertainties. From 
Figure ^, other similar plots for other values of the soft- 
ening and of the initial velocity (not presented here) , and 
the values given in Table we can reach a number of con- 
clusions. 

The dependence of Tioiax on N is reasonably well rep- 
resented by a straight line in the log-log plane and that for 
all values of the softening and /v we tried and for all three 
models. The relaxation time for model D is always smaller 
than that for models P and H. The difference is much more 
important (^^ 1 dex) for the smaller value of the softening, 
than for the larger one 0.15 dex). Similarly, the relax- 
ation time for Plummer distributions is somewhat smaller 
(~ 0.1 dex) than that of the homogeneous one for small 
values of the softening, while for the larger value the two 
relaxation times do not differ significantly. 

Figure || and Table [| also show a trend for the slopes of 
the straight lines. For the small softening the homogeneous 
model has a relaxation time which is roughly proportional 
to the number of field particles in the configuration, and 
that is in good agreement with equation (Q). This is not 
true any more for the more centrally concentrated mass 
distributions, which have a value of Bi less than unity. 
For model P the value of Bi is slightly less than unity, 
but for model D it is considerably so. Thus when we in- 
crease the number of particles in the configuration, we in- 
crease the relaxation time relatively less in more centrally 
concentrated configurations than in less centrally concen- 
trated ones. In other words not only is the relaxation time 
smaller in the more centrally concentrated configurations, 
but also it takes more particles to increase it by a given 
amount. 

These trends for the slopes of the straight lines are 
also present for the larger value of the softening. The dif- 
ferences, however, between the slopes, i.e. between the 
corresponding values of Bi, are much smaller. Thus for 
/v = 1.5 and e = 0.01 there is a difference of roughly 
25% between the slopes corresponding to models H and 
D, while this value is reduced to roughly 8% when e = 0.5. 

Figure || also compares the prediction of equation 
with the results of our numerical calculations for the ho- 
mogeneous sphere. The agreement is fairly good, partic- 
ularly for the larger softening, where the difference is of 
the order of 0.08 dex. It should be noted that these pre- 
dictions were obtained with femax = R, the cutoff radius 
of the system. A yet better agreement would have been 
possible if one used 6max = fR, where /, a constant larger 
than 1. Since, however, this constant is a function of the 
softening used and perhaps also of the value of /v , it is not 
very useful to determine its numerical value. The results 
obtained by using 6niax — I, where I is the mean inter- 
particle distance, are given by open squares. We note that 



Table 2. Coefficients of linear regression in log-log scale 
of Ti-oiax as a function of Vcff 



Model 


e 


A2 


B2 




0-B2 


H 




4.84 


2.72 


0.02 


0.04 


P 


0.01 


4.58 


2.39 


0.02 


0.04 


D 




2.33 


4.28 


0.20 


0.35 


H 




5.20 


2.83 


0.01 


0.02 


P 


0.5 


5.02 


2.53 


0.02 


0.04 


D 




4.42 


2.90 


0.06 


0.11 



they give a bad approximation of the numerical results, 
particularly so for a large number of particles and for the 
value of the softening which is nearest to optimal. For the 
smaller softening the approximation with 6,nax = I is still 
considerably worse than that obtained with fomax = R, 
but the difference is smaller than in the case of the opti- 
mal softening. Whether this will be reversed for an even 
smaller value of the softening or not is not possible to 
predict from the above calculations. Nevertheless, if it did 
happen, it would be for a value of the softening that gave a 
very bad representation of the forces within model H. Thus 
we can conclude that, at least for coUisionless simulations 
which have a reasonable softening, the simple analytical 
estimates presented in section ^ give a reasonable approx- 
imation of the relaxation time if we use 6max — R, but not 
if we use 6max — I ■ The latter gives too high a value of the 
relaxation time, and is therefore falsely reassuring. 

We also compared our results with theoretical esti- 
mates using In A = In(giV) for the Coulomb logarithm 
and d ifferent values of g, as tabulated by Giersz & Heggie 
( 1994|) . We find that they always fare less well than equa- 
tion ([^) with 6max = R, particularly for e = 0.01. Amongst 
the values tried, g = 0.11, proposed by Giersz & Heggie 



( 1994 ), g ave th e best fit for e = 0.5, while the value g = 0.4 
(Spitzer 1987 ) was best for e = 0.01. The differences be- 
tween the results for various values of g are nevertheless 
small. 

To summarise this section we can say that more cen- 
trally concentrated distributions have smaller relaxation 
times and that the difference is more important for smaller 
values of the softening. This argues that the relaxation 
time is more influenced by the maximum rather than by 
the average density. 



5.2. The effect of velocity 

In order to test the effect of velocity on the relaxation 
time we have launched test particles with different initial 
velocities. 

Figure ^ shows the relaxation time as a function of the 
effective velocity of the particles, defined in section U, for 
the three density distributions under consideration. The 
calculations have been made with 64 000 field particles 
and a softening of 0.01 for the upper panel and 0.5 for the 
lower one. The dependence is linear in the log-log plane for 
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Fig. 4. Relaxation time as a function of the effective veloc- 
ity of the test particles, for A^=64 000 and three different 
mass models, namely model H (stars), model P (x) and 
model D (crosses). The straight lines are the correspond- 
ing linear least square fits. The upper panel corresponds 
to e=0.01 and the lower one to e=0.5 

large values of the effective velocity and deviates strongly 
from it for small values. We thus fitted a straight line in 
the log-log plane to the higher velocities estimating for 
each of the mass models separately the number of points 
that could be reasonably fitted by a straight line. We give 
the constants of the regression 

logio(^rclax) = ^2 + -^2 logio(Wcff), 

together with the corresponding error estimates, in 
Table ^ 

In all cases the relaxation time increases with the ini- 
tial velocity of the particles. In order to compare this with 
the analytical predictions of section |^ we note that the 
deviation of a particle from its trajectory due to an en- 
counter should be smaller for larger impact velocities (cf. 
equation |^) , while larger impact velocities imply smaller 
crossing times. Thus from equations (^ and (|^) we note 
that Tieiax should be proportional to the third power of 
the velocity, i.e. that the coefficient B2 in Table ^ should 
be roughly equal to 3. We note that in the homogeneous 
case, which should be nearer to the analytical result, there 
is a difference of less than 10%, presumably due to the fact 



xx 



10" 
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0.01 0.1 
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Fig. 5. Relaxation time as a function of the softening, 
for iV=64 000, /v=1.5 and three different mass mod- 
els, namely model H (stars), model P (x) and model D 
(crosses). 



that the approximations of the analytical approach are too 
harsh. The differences with the results of models P and D 
are on average larger, but strictly speaking, equations ^ 
and (|^) do not apply to them. 

5.3. The effect of softening 

Figure |^ shows the relaxation time Tloiax as a function of 
the softening e for the case of = 64,000 and /v=1.5. We 
note that the relaxation time increases with softening as 
expected. There is no simple linear dependence between 
the two plotted quantities. In fact the relaxation time in- 
creases much faster with log e for large values of the soft- 
ening than for small ones. The point at which this change 
of slope occurs is roughly the same for models H and P, 
and much smaller for model D. In fact in all cases it is 
roughly at the position of the corresponding optimal soft- 
ening, which is roughly the same for models H and P and 
considerably smaller for model D (cf. section O). Thus 
the change of slope must correspond to a change between 
a noise dominated regime and a bias dominated one. 

In the noise dominated regime the sequence between 
the three models is the same as in previous cases. Namely 
it is model H that has the largest relaxation time, followed 
very closely by model P, and less closely by model D. It 
is interesting to note that for high values of the softening 
the three curves intersect. Such values, however, are too 
dominated by bias to be relevant to A^-body simulations. 

Figure ^ shows the only examples in this paper in 
which the error bars are large enough to be plotted, i.e. 
for which CTroiax /^roiax > 0.05. These occur for the small- 
est values of the softening, used here more for reasons of 
completeness than for their practical significance. 
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Table 3. Relaxation time T^oiax of model D with /v ~ 
1.5 and e = 0.01, obtained by the GRAPE-3 (G3) and 
GRAPE-4 (G4) machines. 



N 


G3 


G4 


||G3-G4||/G4 


1000 


1040 


1007 


0.0328 


2 000 


1626 


1633 


0.0043 


4 000 


2 698 


2 707 


0.0033 


8 000 


4 834 


4 845 


0.0023 


16 000 


8143 


8123 


0.0025 


32 000 


14490 


14 548 


0.0040 


64000 


25 858 


25 860 


0.0001 



5.4. GRAPES results compared to GRAPE-4 results 

All results presented so far were made on one of the 
Marseille GRAPE-3 systems (Athanassoula et al. 



199^ . 



Such systems, however, are known to have limited accu- 
racy, since they use 14 bits for the masses, 20 bits for the 
positions and 56 bits for the forces. One could thus worry 
whether this would introduce any extra graininess, which 
in turn would decrease the relaxation time. 

In order to test this we repeated on GRAPE-4, which 
is a high accuracy machine, some of the calculations made 
also with GRAPE-3. For this purpose, we ran 7 configu- 
rations of model D with /v = 1.5, e = 0.01 and different 
number of field particles. The calculated relaxation time, 
Treiax, obtained by the GRAPE-3 and GRAPE-4 runs, to- 
gether with their differences, are shown in Table ||. As we 
can see the results have, in all but one case, a difference 
less than 0.5%. Only in the case of = 1 000 does the dif- 
ference rise to 3%, but as we mentioned before, this very 
low number of particles is hardly used nowadays, even in 
direct summation simulations on a general purpose com- 
puter. 

5.5. Predicting the value ofTj.c\a,x for an N-body 
simulation 

In the standard version of the angular deflection method 
we have used so far all the test particles start from the 
same radius with the same initial velocity. On the other 
hand in any TV-body realisation of a given model the par- 
ticles have different apocenters. It is thus necessary to ex- 
tend this method somewhat in order to obtain an estimate 
of how long a given A'^-body simulation will remain unaf- 
fected by two-body relaxation. Let us consider a simple 
model consisting of a Plummer sphere of total mass 20 
and scale length 9.2, truncated at a radius equal to 30.125, 
i.e. at a radius containing roughly 7/8 of the total mass. 
As an example we wish to estimate the relaxation time 
of a 74668-body|^ realisation of this model which will be 
evolved with a softening of 0.5, a value near the optimal 
softening for model P. For this we will somewhat modify 



the standard angular deflection method in order to con- 
sider several groups, starting each at a given radius. We 
first calculate the relaxation time from each group sepa- 
rately. The relaxation time of the model will be a weighted 
average of the relaxation times of the individual groups. 
The weights have to be calculated in such a way that the 
mean velocities with which the test particles encounter 
the field particles at any given radius represent fairly well 
the encounter velocities between any two particles in the 
A^-body realisation, which is not far from the dispersion 
of velocities. We found we could achieve this reasonably 
well by considering 18 groups, of 10 000 test particles each, 
with apocenters Rmax{i) = 1.25j -I- 2.5, i = 1, ..18. Each 
group starts from a distance such that /v = 0.2 and we 
weight the results of each by 3^*, « = 1, ..18. These weight- 
ing factors were just found empirically after a few trials 
and deemed adequate since they give an approximation 
of the velocity dispersion of the Plummer sphere of the 
order of 10 per cent. It would of course be possible to 
get a better representation by using a larger number of 
groups and e.g. some linear programming technique, but 
since we only need to have a rough approximation of the 
encounter velocities and the fit we obtained is adequate, 
wc did not deem it necessary to complicate the problem 
unnecessarily. 

As expected, we find that the relaxation time and the 
transit time are larger for groups with larger initial radius. 
The range of values we find is rather large. Thus for the 
innermost group we find a relaxation time of 1.6 xlO* 
and a transit time of 4, while for the outermost group the 
corresponding values are 3.6 xlO^ and 51. The weighted 
average of the relaxation time, taken over all orbits in the 
representation, is 3.4 xlO'*, and that of the transit time 
4.9, i.e. nearer to those of the inner groups due to their 
larger weights. 

A comparison with the estimates of the simple analyt- 
ical formula given in equation (7) is not straightforward, 
since one has to adopt a characteristic radius, and the re- 
sult is heavily dependent on that. Thus if we adopt an 
outer radius, where the virial velocity is small and there- 
fore the crossing time large, we obtain very large values of 
the relaxation time, like those we find for the outer parts of 
our model Plummer sphere. On the other hand if we take 
the half mass radius then we obtain Tj-dax = 3.1 x 10^, 
in good agreement with our estimate obtained from the 
weighted average of all groups. 

Our model P is sufhciently similar to the Wc = 5 King 



^ This value of N was chosen in order to have 64 000 field 
particles within a radius of 20 



model used by Huang, Carlberg & Dubinski (1993) to al- 
low comparisons. These authors obtain the relaxation time 
by monitoring the change of energy of individual particles 
in a simulation. They consider only particles which at the 
end of the simulation are near the half mass radius and 
find a relaxation time which, rescaled to our units, is 2.5 
xlO^. Thus this estimate is based on a group of particles 
which have their apocenters at or beyond the half mass 
radius. Applying our own method only to such particles 
also we obtain a relaxation time of 2.3 xlO^, whic h is i n 
excellent agreement with the value of Huang et al. ( |1993|) . 
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6. Summary and discussion 

In this paper we have calculated two-body relaxation 
times for different mass distributions, number of particles, 
softenings and particle velocities. For this we launched test 
particles in a configuration of rigid field particles and mea- 
sured the relaxation time from the deflection angles (mea- 
sured from the theoretical trajectory of the same particles) 
and the transit times. 

We first determine the range of softening values for 
which the error in the force calculation is dominated by 
noise, rather than by bias. These extend to larger values of 
the softening for smaller number of particles and for less 
centrally concentrated configurations, in good agreement 
with what was found by AFLB. We also find them to be 



somewhat larger for models with a larger cut-off radius. 

We confirm that the relaxation time increases with the 
number of particles. Indeed a larger number of particles 
entails a lesser graininess and thus a smaller effect of two- 
body encounters. In particular for homogeneous density 
distributions we confirm the analytical result that the re- 
laxation time is proportional to the number of particles. 
We find, however, that this dependence does not hold for 
all mass distributions. 

We find that the relaxation time depends strongly on 
how centrally concentrated the mass distribution is, in the 
sense that more centrally concentrated configurations have 
considerably shorter relaxation times. This can be under- 
stood if we make an TV-body realisation not by distribut- 
ing particles of equal masses in such a way as to follow the 
density, but by distributing the particles homogeneously 
in space and attributing to each one of them an appro- 
priate mass. Our effect then simply follows from the fact 
that the deviation of a particle trajectory by a single very 
massive particle is larger than that due to a sum of low 
mass particles of the same total mass. 

This dependence of relaxation time on central concen- 
tration is strong. E.g. for a softening e = 0.01 and a /v of 
1.5 the relaxation times of models D and H (or P) differ 
by roughly an order of magnitude. In order to achieve this 
difference by changing the number of particles one has to 
increase them by a factor of 10 also. I.e. in order to avoid 
two-body relaxation ten times more particles are neces- 
sary for a simulation of model D than for one of model 
H, provided one uses the same softening. The difference is 
even larger if the softening is chosen to be optimal in each 
configuration, since the optimal softenings differ by more 
than an order of magnitude, so an extra factor of two is 
introduced to the necessary particle number. 

Also the dependence of the relaxation time on the num- 
ber of particles changes with the central concentration of 
the configuration. We find a shallower dependence for our 
more concentrated models, the difference being more im- 
portant for smaller values of the softening. For a softening 
of 0.01 the difference in the exponent of the power law 
dependence is of the order of 20%. 

We find that the relaxation time increases with ve- 
locity, as expected. The reason is that the defiections in 



two-body encounters are larger when the relative velocity 
is smaller. The dependence of the relaxation time on the 
effective velocity is linear in a log-log plane for the larger 
values of the effective velocity we have considered and de- 
viates strongly from linear for smaller velocities. In the 
simple analytical estimates of section ^ Ticiax is propor- 
tional to the third power of the velocity. Our more precise 
numerical estimates argue that these estimates are only 
about 10 % off for the case of the homogeneous sphere, to 
which they apply. 

The strong decrease of relaxation time with encounter 
velocities entails that two-body relaxation has little effect 
in simulations of "violent" events as collapses or merg- 
ings. On the other hand it may, depending on the config- 
uration, the number of particles and the softening, play 
a role in simulations of quasi-equilibrium configurations. 
Furthermore two-body relaxation will be less of a menace 
in simulations of objects with high velocity dispersions, 
like giant ellipticals which are largely pressure supported, 
than in cases with less pressure and more rotational sup- 
port, like small ellipticals or discs, putting aside of course 
the effects of shape and rotation, which we have not ad- 
dressed here. 

We have also examined the dependence of the relax- 
ation time on softening. We find that, as expected, the 
relaxation time increases with softening. The dependence, 
however, is complicated, and not given by a simple math- 
ematical formula. Nevertheless for the not too centrally 
concentrated models the increase is not too large. Thus 
for our models H and P we increase the relaxation time 
by a factor of the order of 2 if we increase the softening 



by a factor of 10. In this we agree well with Theis (1998). 
The only case where the increase is more pronounced is 
for model D and particularly for high values of the soften- 
ing. It should, however, be remembered that this is a bias 
dominated regime. We can thus conclude that in the noise 
dominated regime the increase of the relaxation time with 
the softening is relatively small. 

Finally we compared results obtained using GRAPE-3 
with those found by GRAPE-4, and found they are sim- 
ilar. From the above results we can deduce that the re- 
laxation times of the two types of GRAPE systems are 
essentially the same. This can be understood as follows. 
Athanassoula et al. ( 199^ ) argued that the limited pre- 
cision of GRAPE-3 does not influence the accuracy with 
which the force is calculated since the error in the calcu- 
lation of the pairwise forces can be considered as random 
(cf. their Figure 5). This is in good agreement with what 
was initially argued by Hernquist, Hut & Makino (1993) 
and Makino (1994), who pointed out that the two-body 
relaxation dominates the error and that the effect of the 
error in the force calculation is practically negligible, pro- 
vided this error is random. 

Our results for the relaxation time are always smaller 
than those given by the simple analytical formula. The 
differences are relatively small if one uses 6max = i? in the 
formula, but quite large if one uses 5max = I- Thus our 
results argue strongly that the former is the right value to 
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use, at least for coUisionless simulations. In this we agree 
with Spi tzer fc Hart ( |197lD, Fa rouki & Salpeter ( |1982D , 
Spitzer ( |1987| ) and Theis ( |199^ ). It should also be noted 
that the analytical estimate obtained with 6niax = I is 
always considerably larger that the numerical result, and 
thus is falsely reassuring. 

By extending somewhat the standard method based 
on the angular deflections we obtained an estimate for the 
relaxation time of an A^-body simulation of a Plummer 
sphere. Comparing it with the results found by Huang et 
al. p9l ) we find excellent agreement. This is ve ry in- 
teresting since the method used by Huang et al. ( 1993 ) 
obtain their estimate of the relaxation time directly from 
an iV-body simulation, i.e. include collective effects. This 
agreement could argue that such effects are not very large, 
and thus gives more weight to the results obtained with 
our simple and straightforward method. 

It is often stressed that galaxies have so many stars 
that their relaxation times are far longer than the age of 
the universe, and thus that TV-body simulations extending 
to long periods of time should have a very large number of 
particles to also achieve sufficiently large relaxation times. 
As a counter-argument one could say that real galaxies are 
not only composed of individual stars, but also of star clus- 
ters and gaseous clouds, which, being considerably more 
massive than individual stars, will introduce two-body re- 
laxation and change the dynamics from that of a coUi- 
sionless system. This, however, is no argument in favour 
of iV-body simulations with short relaxation times, since 
the deviations from the evolution of a smooth stellar fluid 
brought by the graininess of the A^-body system need not 
be the same as those brought by the compact objects, star 
clusters or gaseous clouds. It is thus necessary in A^-body 
simulations to strive for high relaxation times and believe 
the results only for times considerably shorter than that. 
If desired, the effects of the compact objects, star clusters 
or gaseous clouds can then be studied separately. 
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